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Abstract 

A new analytic expression is derived for the resonant drive of high n Alfvenic modes by particles 
accelerated to high energy by Ion Cyclotron Resonance Heating. This derivation includes finite 
orbit effects, and the formalism is completely non-perturbative. The high-n limit is used to calculate 
the complex particle response integrals along the orbits explicitly. This new theory is applied to 
downward sweeping Alfven Cascade quasimodes completing the theory of these modes, and making 
testable predictions. These predictions are found to be consistent with experiments carried out on 
the Joint European Torus [P. H. Rebut and B. E. Keen, Fusion Technol. 11, 13 (1987)]. 
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I. INTRODUCTION 



The behaviour of highly energetic ions in tokamak plasmas is one of the most important 
issues in the physics of burning plasmas^, and thus of great interest for both theoretical 
and experimental study. Alfvenic instabilities which can tap the free energy in the energetic 
ion distribution^ are often excited in present-day tokamak experiments. This is due to the 
resonant interaction between shear Alfven waves and ions either accelerated by ion-cyclotron 
resonance heating (ICRH)31 or produced by neutral beam injection (NBI)R It is anticipated 
that similar Alfvenic instabilities will be excited in burning plasmas such as ITER by the 
fusion-produced alpha particles 1 2 . The destabilization of Alfven Eigenmodes (AEs) typically 
results in spatial redistribution (potentially global) of the fast particles^, which may lead to 
quenching the plasma burn and damaging the first wall. More recently AEs have also been 
used for diagnostic purposes^, divulging information about the plasma interior inaccessible 
to other diagnostics. 

In order to determine the potential effects of AEs and what information they can reveal 
we require a complete model for both the modes and their interaction with the energetic ions. 
In the analytic study of energetic particle driven instabilities the focus has often been on 
the excitation of bulk plasma eigenmoded 7 1 8 1 9 1 10 -! There is also a lot of interest in instabilities 
associated with non-perturbative energetic particle modes^^. We continue this progress 
here by developing a new analytic model for the resonant drive of shear Alfven perturbations 
due to ICRH-accelerated ions, suitable for analyzing both classes of instabilities. 

The distribution of ions created by ICRH-acceleration is a distribution of particles with 
predominantly perpendicular energy and sharply peaked in pitch angled The spatial dis- 
tribution of these particles is such that all the tips of the banana orbits lie on the surface 
where the local ion cyclotron frequency resonates with the applied RF wave. In a typical 
tokamak this is an approximately vertical surface, R ~ Rrf, and we classify the heating 
by where Rrf lies compared to the magnetic axis R . If Rrf < R we have high-field- side 
heating, Rrf > Ro gives low-field-side heating and Rrf = Ro is known as on-axis heating. 
In a typical discharge with ICRH the flux-surface-averaged density of fast particles can be 
hollow for low-field-side heating, and is usually centrally peaked when there is high-field- 
side or on-axis heating. Given the energetic particle distribution we can calculate the fast 
particle effect on any given mode. 
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Our considerations will be limited to the case of Alfvenically-polarized fluctuations with 
high toroidal mode number (n ^> 1). In Sec. [n]we use this approximation to amalgamate 
previous results and produce a linearized MHD description of the problem including a kinetic 
contribution from the energetic particles through their perturbed anisotropic pressure. 

In Sec. Ill , buliding on previous work^^HS ) we cas ^ wave-particle resonance in 
Hamiltonian form and explicitly determine the transformation to a set of action-angle coor- 



dinates for the complex bounce-precessional energetic particle orbits. In Sec. [TVJwe use these 
coordinates to express the linearized drift kinetic equation for the fast particles in a simple 
form. We then solve this equation for the fast particle response and use it to calculate the 
contribution that the fast particle resonance provides to the Alfven wave equations devel- 
oped in Sec. |TT} The result is expressed in terms of complex integrals over the unperturbed 
orbits, to make further progress we focus on the case of radially extended perturbations, for 
which 

dm , 

7T < - 1 
or r 

(where r is the minor radius of the plasma and m is the poloidal mode number). This 
approximation is valid for high-n Alfvenic perturbations if the magnetic shear is sufficiently 
low over a substantial fraction of the minor radius. The radial structure of the mode is 
known to be sensitive to the radial variation of the g-profile. The Alfven frequency variation 
is itself controlled by the g-profile variation which is amplified by the large mode number. 
Thus if the low shear region extends for multiple poloidal wavelengths the relation in Eq. 
([T]) should hold. The envisioned mode structure is shown schematically in Figure. [I] 

This assumption and the thin orbit approximation imply the mode is much broader than 
the orbit width, which alllows us to average over a radial distance wider than an orbit width 
yet smaller than the mode width. The sheared magnetic field configuration leads to the 
result that the energetic particles experience different phases of the mode at different radii 
along their orbit; as these differences are multiplied by the large mode number n they lead 
to highly oscillatory contributions to the integrals along the energetic particle orbits. The 
averaging procedure described allows us to calculate the integrals along orbits in closed form 
by only considering the non-oscillatory contributions in the vicinity of the turning points. 
This enables us to formulate a local expression for the resonant particle response, which can 
then be used in a variety of problems associated with ICRH-driven Alfvenic modes. 

Many modes satisfy the condition in Eq. 0, including Low-Shear Toroidal Alfven Eigen- 



FIG. 1: Schematic mode structure, depicted by one equipotential surface, superposed with the 
poloidal cross section of an energetic particle banana orbit 



modes^ 17 * 18 *, which appear as "tornado" modes in experiment ^ 19 * 20 ! Eq. ([!]) is also applicable 
to the Alfven Eigenmodes observed in reversed shear discharges on JET^H, JT-60U^ and 
Alcator C-MocP^. In order to demonstrate the utility of our technique for performing res- 



onant drive calculations we apply it in Sec. V to the problem of the Alfven Cascade (AC) 
modes in scenarios with weakly-reversed shear. An explanation for the downward sweep- 
ing mode observed in such scenarios as a weakly-damped propagating quasimode has been 
proposed^, which we complete by calculating the radiative damping and including the res- 
onant contribution from the energetic ions. The damping and the drive can be evaluated 



for both downward and upward sweeping modes. We culminate in Sec. |VII| by comparing 
our theoretical results with experimental observations from the JET experiment!^. 



II. ALFVENIC PERTURBATIONS IN A TOKAMAK 

Motivated by fusion applications, we will make the following assumptions from the outset. 
Firstly, in all cases /3 (the ratio of thermal pressure to magnetic pressure) will be considered 
to be vanishingly small (/3 <C secondly we look at modes elongated along the field line, 
i.e. having parallel wavenumbers k\\ = (qR)^ 1 (m — nq) much smaller than perpendicular 
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wavenumbers k± « m/r, finally we will make the approximation of small inverse aspect 
ratio (e = r/i? < 1) and shifted circular flux surfaces where necessary. 

Previous works on AEs have computed the contribution to the linear MHD equations 
governing Alfvenic perturbations due to an anisotropic pressure^^, due to pressure gra- 
dients in the thermal plasmaP^ and due to coupling of the shear Alfven perturbations to 
acoustic and compressional Alfven perturbations^. Assuming a background plasma with 
magnetic field B, pressure P and density p and that the toroidal mode number n of the 
mode is large we can combine these results to give the following wave equation: 



1 - B ■ V I tLv^ B ■ V$ 



v\ x dt 2 V B 2 



Ze 
c 



(2) 

d 3 vv D -V5f + -^bx k- V[(6x VlnP-7^Vx b) ■ V$] , 



where va = B/\Z4np is the local Alfven speed, b = B/B the magnetic field direction, P the 
thermal plasma pressure, n = b ■ Vfo the curvature, Vp = jj^ (uVB + Mv 2 b x b ■ X7bJ the 
fast particle drift velocity, \x = Mv 2 _/2B the magnetic moment and 7# the ratio of specific 
heats of the background plasma. The fast particles' perturbed distribution function is 5f 
and their mass is M, charge Ze and cyclotron frequency Q = ZeB/Mc. We have followed 
Ref. [29] and expressed the electromagnetic perturbation as 

^ = Vx((bx V$) x b) , (3) 
1 9$ 

SB =-^w < 4 > 

This formulation of the problem relies on (3 1 and k\\ k± to eliminate the coupling to 
the compressional Alfven wave (fast magnetosonic mode), and on 2q 2 ^> 1 to neglect the 
coupling to acoustic perturbations (this is equivalent to assuming that the Geodesic Acoustic 
Mode (GAM) is supersonic). This decoupling is carried out in detail in Ref. 

To compute the fast particle resonant contribution to Eq. ([2]), we express the integral over 
velocity space as a curvature coupling between the Alfvenic perturbation and the perturbed 
anisotropic fast particle pressure (again using the fact that f3 is small to simplify the drift 
velocity), 



Ze f , Q _ „ „ 47r 
c / " ' ' B 2 



d 3 vv D .V6f=—bxK- WSP ± , (5) 
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where we have neglected the small parallel pressure in comparison to the perpendicular 
pressure. We now merely have to complement Eq. ^ by expressing the perturbed fast 
particle pressure 5P± in terms of $. This has already been done in Ref. for the non- 
resonant component of SP± , thus we only need consider the resonant component here. 



III. HAMILTON! AN FORM OF THE UNPERTURBED ORBITS 

In this section we reduce the dynamics of the particles to a completely integrable Hamil- 
tonian form, keeping finite orbit width effects but neglecting finite Larmor radius (FLR) 
effects. We start from the Hamiltonian form of the Littlejohn Lagrangian 30 as presented in 

Ref.na 

L = p e e + p^ + p c c - h(p v , p e , e, p c ), (6) 

Pe = ^ Ae + Mv ^ (8) 

C D 

P r= ^A v + Mv,^-, (9) 



n-^B + ^U-^A.) (JL). (io) 



Ze If 

<:\l - ' 2M V " c "*J \B V/ 

Here subscripts denote the covariant component of a vector quantity. In addition ( is the 
gyroangle about the equilibrium field and we have picked orthogonal coordinates r, 9, if 
corresponding to flux surface label, poloidal angle and toroidal angle respectively. The 
definition of these coordinates is chosen such that B r vanishes (equivalently r is orthogonal 
to 9). We have also picked a gauge such that the unperturbed A r vanishes. This is critical 
to being able to write the v ■ x term in the Lagrangian in Hamiltonian formP^EH. j s 
understood that all quantities should now be viewed as functions of P^, Pg, P v and (, 9, if as 
one can define r and v\\ implicitly from Eqs. (|9| and (|8j). 

In order to express the motion in action angle form, we follow Ref. HH and introduce 
a generating function g for the transformation to new canonical coordinates <f, 9, £ and 
momenta Pp, P§, P^. We now consider a function g that leaves P v and P^ unchanged i.e. 

g = i> + P~ C C + g*(P?; Pg; Pf, 6), (H) 
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which gives the following equations for ip and 9: 



The function g* is determined from 



In order to solve this we need an explicit form for Pq as a function of 9 and conserved 
quantities P$, P§, P^ or equivalently P v , and P. This nessecitates picking an explicit 
form for both A and B subject to the orthogonality constraint mentioned above. 

In a large aspect ratio tokamak we can do this by expanding all quantities in the inverse 
aspect ratio to fincP} 

r 

/vdv 
WY (15) 



A e = -B Q r 2 -B a(r)cos9, (16) 



B<p = B R , (17) 

2 

'Pog(r) 



fi e = ^o^r fl - ^costf ) , (18) 



with the function a defined by 



r 

.2 



a ( r ) = / I +rA' + r [ ^-ds I dr, (19) 



R 



o 







with primes denoting radial derivatives, A(r) the Shafranov shift , s a dummy variable, Po 

the magnetic field strength on axis and Rq the major radius of the magnetic axis. We also 

note that there is a simple transformation from orthogonal coordinates to straight-field-line 

coordinates r, 9, (p which are convenient for describing Alfvenic modes (see Ref. I3T]) . 

a'(r) 
r 

<p = <p. (21) 
In what follows we replace P^ by a new radial coordinate f such that 

P ' = -T*JWY (22) 
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We thus have Pe given by Eq. (J8J) as a function of r, 9 and v\\ and must use Eqs. (22) and 
Q to eliminate r in favour of f, P^ and 6. We start by neglecting the FLR correction in 
Eq. (|8]), finding that Pe is now a function only of the spatial coordinates, 

Pe « (23) 

We now make use of the fact that we are dealing with ICRH-accelerated particles, which 
are poloidally trapped between bounce points where v\\ changes sign. The majority of these 
particles perform banana orbits where they deviate only slightly from the flux surface of the 
bounce points, thus we write, 

r = f + 5, (24) 

where 5 is the deviation from the flux surface and we identify f as the radius of the bounce 
points. We will now assume r,f ^> 5. This means that in most places we can replace r with 
f, as effects of the order of 5/r are the smallest effects we wish to keep (finite orbit width, 
but not FLR effects). However in doing so we limit ourselves to orbits that do not pass close 
to the magnetic axis as r ~ 5 there, and such orbits are no longer banana orbits but so-called 
potato orbits. The thin orbit approximation allows us to find S explicitly and eliminate r 
in favour of f and S as required. In the potato orbit regime a different approximation must 
be used to eliminate r. As the AC modes under consideration here are well away from the 
magnetic axis we make the thin-orbit approximation in all further analysis. 
We thus obtain, 

_ Roq(f)v\\ 

r = r-\ — (25) 

r sZ 

We also introduce the bounce angle 0&, being 6 at the v\\ = point, which we will use instead 
of Pg. Eliminating r we find, 

P e = -^[ -B f 2 - B a{f) cose) + q(f) y^J ~^ P C B o ^cos^-cos^. (26) 



We can now solve Eq. ( 14 ) for g* by substituting the following ansatz, 



g* = ^ Qf 2 fl - a(f) sin ^ + q(f) ^R^yj^-P^G, (27) 



where G now satisfies 



~d6 



cos 9 - cos 6 b . (28) 



Using Eq. ((271) in Eqs. (Il2|) and (fl3l) results in: 



qa 

m = m - q Q + 1— sinO - G(9)xi(f) + 
r 



9 = Ugrq, 



2M fR 

f 



^ - X2(f) ) 9, 

3/2 



1 dG 
sin 9h d6 h ' 



(29) 
(30) 



where the functions Xi( r ) an d Xz{ T ) are defined as follows, 



Xi 



q(f) d / 2P C 



f dry MVL 

P^ cos 6b 



X'2 



rR q 2 (r) 
q(f) 



Mu s R Q f 



(31) 
(32) 



The general solution of Eq. (28) is, 



G = ±2 3 / 2 



sin 



E 



CSC 



+ A, 



(33) 



with E(z\m) the elliptic integral of the second kind with modulus m 32 , and A an arbitrary 
constant. Choosing our orbit to start from 9 = —6b, 6 = — n, pass through 9 = 6b, 6 = 0, 
and end at 9 = —9b, 9 = n determines the constants and the ambiguous signs, giving: 



G - 
G - 
and 



2 3/2 sin 
2 3/2 sin 



CSC 



CSC 



-E 
-E 



CSC 



CSC 



when 9 G (-vr,0), (34) 
when 9 E (0, ix), (35) 



d9 


R q 

\ v \\ 


dd ~ 


d9 


R q 


d9 ~ 


h 



when 9 G (0,7r), 



with 



'II 



2B f Ze dG 



Pr 



(36) 
(37) 

i? M 2 c^dB' (38) 
This completes the description of the transformation to the coordinates 9, (p from the or- 
thogonal coordinates 9, ip. 

To complete the analysis of the orbits, we must determine the frequencies with which 
the particle performs its periodic motion in 9, (p. These frequencies however are merely the 
particle bounce frequency and toroidal precession frequency which in the large aspect ration 
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limit are given by£ 



Uf ~ 2Roq{f)K{K)y (39) 

_ c Mr) f + (k' - 1) AT(«) ggOO-jTOQ l 

^ ~ Ze i? f \ 4Slrj K(k) + K(k) J l4Uj 

with k = sin(6 l fc/2) and if and £7 the complete elliptic integrals of the first and second kinds 
respectively. We also write the magnetic shear as 

f \ r d<1 (AV\ 

q[r) ar 

IV. CALCULATION OF THE FAST PARTICLE RESONANT RESPONSE 

Having found action angle coordinates in which the unperturbed motion of the energetic 
particles is completely integrable, we now turn to the problem of calculating the energetic 
particle response to shear Alfven perturbations. Similar derivations for perturbative modes 
and deeply trapped particles have been done before in Ref . [TO] and Ref . [TJ)J For the case we 
will consider in Sec. |V]we expect the resonant instability drive to be small compared to the 
wave frequency and thus for simplicity we will only consider the leading order contribution 
to the resonant fast particle distribution function in the current section. We start by writing 
the collisionless kinetic equation for the energetic particle distribution function as 

ft- - < 42 > 

with the time derivative being taken along the particle orbit. This can then be linearized 
by expanding the distribution function as / = F + 5f where 5f is the perturbation due to 
the mode and F is an equilibrium distribution function that only depends on the constants 
of motion. The resulting equation is 

with the derivative on the left evaluated along the unperturbed orbits, and the one on the 



right along the perturbed orbits. We now write the left hand side of Eq. (43) using the 



coordinates found in Sec. Ill and average over the fast motion in £ to obtain 



d Sf=(i + ^ + "el s )Sf- (ID 



dt J \dt *d(p °86 
10 



If we pick F = F(U, fi, P$), where U = \mv 2 is the particle energy, then the right hand side of 
Eq. (43) can be evaluated by conventional means^ 34 * 35 ! for the electromagnetic perturbation 



given by Eqs. ^ and Q, which results in 

d „ Zed (d§\ OF Ze B$dF . . 

d~t = ~~c\p \~dt ) dP^ ~ ~c Vd ' ~didU 1 (45) 

where derivatives of F are taken holding fi and either U or Pg constant. We note that the 



second term on the right hand side of Eq. (45) is smaller than the first term by a factor of 



r/R and thus we neglect it and focus on radial gradients of the energetic particle distribution 
as the main instability driving mechanism. Combining Eqs. (44) and (45) results in 

(I + <4 + "4) V = - Z i^^-% (46) 

where we have made use of a Fourier expansion of $ in time and in straight-field-line coor- 
dinates (0,0) as follows: 

$ = <& m {r)e in *- im °- iu}t . (47) 



Axisymmetry implies (due to Eq. (|29|)) that Fourier harmonics in (p are decoupled, and thus 
one harmonic only gives rise to one (p harmonic. Expanding Sf in a series of 6 harmonics 



and a single if and time harmonic then gives an immediate solution of Eq. (46) : 



. sr^ Ze^ m (f)K(l,m)e in ^- ife - iu;t dF 

of = —inuj y : — . (48) 

c uj — nuj(p + lujfj oP v 



l.rn 



To transform from r, 9, (p straight-field line coordinates to f, 9, <p~ coordinates, we have made 
use of the slow radial variation of $ m to replace $ m (r) by $(f) and defined the internal 
integral K, which transforms straight-field-line Fourier harmonics to action-angle Fourier 
harmonics, by 



2tt 



K(l,m) = ^ J dp J d9e in *- imS e- in v +il§ . (49) 

-7T 

In general, the evaluation of K involves a complex integral over the entire orbit. However, 
in the high-n approximation, with fen <C fej_, we can evaluate K asymptotically, which we do 
in Appendix [A] 



We now use Sf as given by Eq. (48) to compute 5P for inclusion in the mode equation 
(Eq. Q), via Eq. ([5]). In order to integrate 5/, we write integrals over all velocities as 
integrals over U, /i and (. However it is important to note that whilst in the definition of 
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the pressure the velocity integration is done at constant r, 5f is a function of f. We make 
use of the thin orbit apprixmation and slow radial variation of all quantities to replace all 
instances of f with r in the integrand, and perform all integrals at constant r. Thus 



MR q(f)ue 



de 



de 



(50) 



where Eq. (36) has been used to replace |t>ii|, and a 
for 5P\ to 



This simplifies the expression 



5P ± = 2tt J dUd » 



de 



de 



Subsituting 5f from Eq. (48) we then have 

B f 

5P ± = - 2mnuQ ]^^ 1 $m / dUdfi 



Lm.i 



MiR q(f)LUg 

H K(l,m)e in ^- il§ - iuJt dF 
ujq u — nujq, + lujg dP^ 



de 



de 



(51) 



(52) 



We now expand 5P± as a Fourier series in (p and 9, and interchange integration over (p and 
9 with integration over U and fi to find 



SP± 
5P ±j 



re 



B, 



intp—ij9—iujt 



—2itinujVl ° — $ m / dUdu, 
MR q ' I 



fi K(l,m)K ext (l,j) OF 

"tUg U — nUJcp + lujg OP^ 



(53) 
(54) 



l,m tv-,. 

We have defined the external integral K ext , which is similar to K except involved with the 
inverse transformation back to straight-field-line harmonics, by 

2"7T 7T _ 



1 

47T 5 





de 


1 'dtp j ' de 


de 



^intp—iW g—irup+ijB 



(55) 



and again defer its evaluation to Appendix [Aj 

We now simplify SP±j by selecting only the resonant contribution, i.e. we use the formula 



(56) 



dx f f dx 

= -i dx5(x)+V / — , 

X J J X 

with V j dx denoting the Cauchy principal value of the integral. We define 5P± to be the 
part of the perturbed pressure associated with the resonant delta function. We then use 
the nature of the background fast-ion distribution F to perform the integration over U; the 
approximation that all the fast ions bounce at the resonant magnetic field strength (where 
the ICRH deposition is localized) is equivalent to factorising F into 



F = F Q {P ¥ ,^U)8{U-iiB res ). 



(57) 
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We finally combine Eqs. (|54| and (57) into 



8P% = -2Tmutt— ^-Y> m / d^S(u;-n^ + lu; s )-^K(l,m)K ext (l,j)-^, (58) 
MK q ' J LO e Or^ 

with all instances of U evaluated at the resonant energy fiB res . If we average this over a 
small radial region (larger than the orbit width but less than the radial wavelength of the 
mode) or equivalently neglect components of 5P that vary rapidly in radius we can use the 



explicit expression for K(l)K ext (j) from Eq. (A21) to ultimately obtain 

5P R ±J = -2W)— £ $ m / d/iS (w - nw p + lug) ii—^-^., (59) 



if # b < I or 



2 



" ■ cos ((m — jjy ) 



= -2imuVL y2&™ djj5 (to - nu$ + lug) 

MK q ~ J 

l,m 







a;e7r|n(l7*| dP^ 

(60) 

if #b > |. All quantites in these expressions are now functions of r and not f (see discussion 



after Eq. (49)). The quantities G* and 6* are defined by 

2P ? Qfxi 



G* 



M q 2 R 3 ooj 
1 
2 



f- sin fl* v 7 ^^* - cosfl 6 ) , (61) 



cos0* = ^-cos0 6 . (62) 



V. APPLICATION TO ALFVEN CASCADE DRIVE 

In this section, we apply our analysis to the problem of Alfven Cascades in shear-reversed 
discharges driven by ICRH-accelerated ions. We expect the resonant ions to have signifi- 
cantly lower energies than the non-resonant ones. We therefore use the small orbit width 
approach detailed above for the resonant ions and the large orbit width limit considered 
in Ref. [26] for the nonresonant ions. Alfven Cascades can be adequately described by a 



single Fourier component $ m of the perturbed field (see Eq. (47)); taking the corresponding 
component of Eq. ^ results in: 

dr \v\ " v 2 A ) dr r 2 \v\ " v 2 A v\ 



in 



4nem dnoh 2tt / m R R , d , R R A 

P^^T $m + iB^Ra \7 {6P±m+1 + 6P±m ~ l} + dr ( 6P±m+1 ~ 5P±m - l} ) ■ 

(63) 
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We have introduced the geodesic acoustic frequency uJq 



MB 2 



(T e + jTA and the pressure 



gradient induced frequency shift uj\ 



MR 2 ' dr 



r4- (T e + Ti) in accordance with Ref. (24] and 



Ref. 



We can write the resonant pressure contributions using Eq. (58) as: 
1 



-2-KnujQ- 



B 



2MR q 



$ m / dfi5(uj — mom + Iujq)— KM, m) 



dF 

dP ln 



dP Lr , 



2 {^L + i " 5PL-i) 

= -2nnuQ ^° $ m [ dfi5(tJ - nu^ + lug) — K 2 (l, m" '° l " 
2MR q t-f J lo§ 

K 1= J2 K ( l > m)K ext {l, m + s), 

s=±l 

K 2 =J2 sK ( l > rn)K ext {l, m + s). 



(64) 

(65) 

(66) 
(67) 



s=±l 



With these expressions Eq. (63) reduces to a second order differential equation for the shear 
Alfven perturbation $ m , 



d 



2 2 \ 

UJ -u A ) 



uJ 



U A — UJ 



2 \ m 



iD&m + D 2 



dr v AJ dr v ~ A n> r 2 "* ' dr 

where we have abbreviated the fast particle terms by defining D\ and D 2 as follows: 

2nv\m , R R >. — 



iBqRq r 



9 
CP 

or 



(68) 

(69) 
(70) 



iB Q Rodr v" ±m +! 

and defined frequencies of the Alfven continuum (ua) and the offset (uh) from the continuum 
as follows: 



2 ; 2 2 , 2 



a; 



H 



Gi 

u r dn oh 2 

— a hoV 

m n,- ar 



(71) 
(72) 



Before engaging in a detailed analysis of Eq. (68), we can immediatly determine several 



general features of the modes. Making a local approximation, and neglecting the small 
radial derivatives we have 



or 



2 2 \ 

U) A -UJ H ) 



or 



(73) 



where A is a positive constant, and K\ is defined by Eq. (A22). This leads to the conclusion 
that the fast particle contribution is destabilizing if: 
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• 9b > | and ^ < 0, i.e. High field side ICRH and a peaked radial profile, 

• Q\) < | and |£ > 0, i.e. Low field side ICRH and a hollow radial profile. 

As these are the expected profiles of ICRH ions, we expect the mode to be observed indepen- 
dent of heating location. We also note that these results are for n > 0. Looking at the exact 



expression in Eq. (64) we see that the destabilizing term is proportional to n. Therefore 
we do not expect any modes with negative n for normal profiles of ICRH-accelerated ions. 
Returning to Eq. ([68]), we proceed as in Ref. [2U and use our assumption of small radial 



gradients (Eq. ([!])) to determine that the second term on the left hand side of Eq. (68) is 
larger than all the other terms. This results in 

u 2 « u 2 = u\ + u 2 H , (74) 

i.e. the real frequency of the mode must be close to the local Alfven frequency with the offset 
determined by the nonresonant fast particle density gradient and plasma pressure gradient. 
This means that we can neglect the small spatial dependence in the radial derivative term, 
d 2 , m 2 uj 2 - u A (r) - uo 2 H (r) , 1 , . . d , 

or z r A ^hvo) iO H \r§) \ or 

We order the radial dependence of u A (which is predominantly from the radial variation of 
q) through 

t 2 d 2 

y A (r) - u 2 A (r )\ ~ uj H (r ) — — . (76) 

As uju "C uja, this means that we can neglect the spatial dependence of uj 2 h everywhere. 
Expanding about q = qo, the value of q at the reverse shear point, where r = ro , we find 
(writing u 2 . = 0/4 (ro) + U}jj(ro) to the required accuracy) 

' -$ m = (A - V z 2 - z 4 ) $ m + iE$ m + F—$ m , (77) 



/ 2 // 2\ 1/3 

r — r 1 m Va%Tq x 



A 



r V 2u H q£R { 

2 , ,2 / /), ,2 „4 d2^2\ 1/3 



(78) 



, = ( ^- m) UsW ' (80) 

£ =t(*^) 2/3 A, (81) 

1/3 

> > „ , ^2. (82) 



r ( 2q 2 l u H R 
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We then introduce the auxiliary quantity \& = e Fz,/2 $ m , which reduces Eq. (77) to 

d 2 



dz- 



A + iE - -F - nz 
4 



2 _2 .4 I ^ 



(83) 



We also find it convenient to reintroduce explicit time dependence by writing $ r 



&m{t)e luJot with $ m (t) a slowly varying envelope. This results in 



.or 



We can then rescale time by t = t 

. 9 ~ 



2/3 



lu 2 ) * « 2^ |^. 



to obtain the final mode equation, 



, d 2 

df dz 2 



* + (r^ 2 + z 4 ) * 

* = e 



(84) 
(85) 



which is a time-dependent Schrodinger equation with a quartic potential hill for large z and 
a local hill or a well for small z, depending on the sign of 77. The expression for c^o shows 
that a downward sweeping Alfven Cascade will have rj > and an upward sweeping cascade 
will have rj < 0. Thus we naturally see that the upward sweeping modes are generally real 
eigenmodes localized in a potential well and the downward sweeping modes are decaying 



quasimode solutions of Eq. (84) atop a potential hill. The transformation to \I> means that 



if we find a solution of Eq. (84) with damping rate 7(77) we have a solution of Eq. (75) with 
damping rate T given by, 



D x 1 / uJ 2 H v A q'o 



"„2 \ 2 / 3 



+ 



m 



2 2u; 2ujq \ 2mq 2 R 



Tfa) 



(86) 



VI. ANALYSIS OF EQ. (84) 



Having reduced the problem to the simple form of Eq. (84) we now proceed to analyse 



this equation to find 7(77). Our goal is to find the least damped solution of Eq. (84) in the 
absence of a drive. This corresponds to the most unstable/least-damped solution for $ m 



through Eq. (86) in the presence of drive. Examining the form of Eq. (84) we see that for 



large \z\ the solutions are asymptotically WKB solutions of the form, 



~ — e 

z 



±iz 3 /3 



(87) 



which are radially ingoing / outgoing waves. By analogy with the Schrodinger equation the 
z-flux of |^ 2 | is given by i (^fd z ^f* — ^/*d z ^, which shows that the outward flux of |^ 2 | 
at large z is proportional to ^ 2 |\I /2 |. We therefore expect perturbations near the origin to 
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last longest before radiating away to large \z\. It is also seen immediatly from the WKB 
solutions that as the energy is radiated away from the origin, it quickly reaches ever shorter 
scales and thus formally our assumption in Eq. Q becomes invalid. However, if we assume 
that there is some (presumably kinetic) dissipation at small scales then we can envisage a 
situation in which the region around z = is well described by our model, and it radiates 
energy away at a rate that is independent of the exact small scale dissipation mechanism. 
In order to calculate the characteristic lifetime of a perturbation we use an initial value 



code to solve Eq. (84). We discretise both in space and in time, using a standard Crank- 
Nicolson scheme. In order to mimic a generic small scale damping mechanism at large z 
at each timestep we multiply the solution by a smooth windowing function that artificially 
damps fluctuations far from the origin. This provides an absorbing boundary away from the 
origin allowing energy to radiate away from the origin realistically. 

As this is a Schrodinger equation with a potential hill, we expect there to be a complete 
set of solutions with outgoing wave boundary conditions and quantized damping rate. Thus 
any initial perturbation can be expanded in these functions, and after a few damping times 
only the least-damped of these solutions (quasimodes) remains. The resulting damping rate 
is thus independent of initial conditions. We find this damping rate by taking the absolute 
value of ^ at the origin and fitting ln(|\l/|) to a linear function of time, after removing the 
initial transient phase where many modes are still present. 

We also seek an analytic approximation to this damping rate in order to validate our 
code. For large positive rj the problem can be transformed into the standard anharmonic 
oscillator with a complex quartic term. This problem has been studied before, and so we can 
use standard results^ to obtain for large rj 

1{V) = -Vv(l- l*V~ 3/2 + §?T 3 ) • (88) 
However, as we only require the real part of 7 we can neglect the ^~ 3 / 2 term. We note that 



the analysis in Ref. [2U of Eq. (84) for large 77 contains misprints and should read (in our 



notation) that \I/ ~ e l ^ f,z I 2 with A ~ i-Jrj, in agreement with Eq. (88). As a further check we 
confirm that in the case of 77 = we reproduce the ground state of the anharmonic oscillator 
correctly^. In Fig. [2] we present the analytic approximation for 7 with the numerical results 
obtained by our finite difference solver. We note that the exponential tail for negative rj 
gives the small radiative damping of the usual upward sweeping eigenmode. 
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FIG. 2: Numerical solution for 7(7?) compared to analytic result for large r/ (color online) 

It is important to compare our results to alternative interpretations of downward sweeping 
ACs. Whilst we consider a weakly-radiatively-damped mode atop a potential hill, it is 
possible to construct solutions where kinetic effects provide a potential well for standing 
wave ^38139] v j a coupling to kinetic Alfven waves. 

We work in the same regime as Ref. EH] where T, given as the sum of the drive and the 



continuum damping by Eq. (86) is not vanishingly small. In this regime, the discrete kinetic 
modes have a spectrum consisting of multiple eigenmodes with a frequency spacing 5u given 
byP 




We compare this to the damping rate 7^ of our longest lived mode in the limit 77 3> 1, i.e. 
away from the lowest frequency of the cascade or for very small q", 



1± 




Thus, if the condition 



Pitn uj a 



-1/2 

1 ) . (90) 



UH > 1 (91) 
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is satisfied, we expect the resulting behaviour to be well described by our continuum approx- 
imation. This is because the radiative damping is in reality due to phase mixing that arises 
from the many different kinetic modes that make up the initial perturbation. The discrete 
nature of these modes is only apparent on timescales longer than Suj^ 1 by which time the 
mode has either decayed to zero or grown beyond the applicability of linear theory. In the 



opposite limit, where Eq. (91) is violated, the results of Ref. [39] have been extended in Ref. 



to find a discrete kinetic mode with a broad radial scale. 



We now consider the condition given by Eq. (91) for typical experimental parameters. 



Using Eq. (72) and neglecting the pressure gradient allows us to estimate uj h by 



^ - \/--^f , (92) 
uj a V m uj A rii L h 

where Lh is the density gradient scale length of the fast particles. For typical JET parameters 
fSL ~ io 3 , m « 10, ^ « lO" 3 , a « 100 and f- « 3) we find 



and 



UJa 2 



r uj h 



5 > 1, (94) 
Pirn uj a 

so that the continuum-like description for the mode is suitable for the problem under con- 
sideration. Near marginal stability, i.e. for the case 7l — 7d ^ Id, where jl is the linear 
drive due to resonant particles, an extra factor lL ~^ d appears in Eq.(94) , which may violate 
the applicability condition for the quasimode approach and a more delicate description of 
the discrete spectrum would be required. 



VII. COMPARISON WITH EXPERIMENT 



In relating our results to experimental data, we need to highlight the fact that the down- 
ward sweeping ACs are very rare compared to upward sweeping ACs. Since 1997, when 
magnetic spectrograms became available for every JET discharge, downward sweeping ACs 
have been detected in only 26 discharges as opposed to more than 5000 discharges where 
upward sweeping ACs have been observed ( magnetic probes, O-Mode inteferometryS^ and 
X-Mode reflectometrySll contribute to this data). Due to the scarcity of downward sweeping 
ACs, some essential plasma diagnostics were unavailable at the time of their observation. 
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In particular, MSE measurements^ of the q(r)-profile were only made in 5 discharges out 
of the 26. Although 26 examples of downward sweeping ACs provide very limited statistics 
they are still indicative of why downward sweeping ACs are so rare. 

First, there is no obvious correlation between downward sweeping ACs and the type of 
energetic particle population in the discharge. The downward sweeping ACs were observed 
in discharges with NBI only (3 discharges), combined ICRH and NBI (9 discharges), and 
ICRH only (14 discharges). In the cases of plasmas with ICRH, high-field side (up to 
Rrf — Rq ~ —30cm in discharge #56947, t = 3.8 s), low-field side (up to Rrf — R® ~ 35cm 
in discharge #65550, t = 5.2 s), and on-axis ICRH were used. It is therefore unlikely that 
downward sweeping ACs require an extraordinary energetic particle distribution. 

Second, downward sweeping ACs are seen in discharges with magnetic fields ranging from 
1.8 T (discharge #72691) to 3.45 T (discharge #46863). Although no measurements of ion 
temperature exist for most of the downward sweeping ACs, the variation of the magnetic 
field, together with the difference in the types of heating, is likely to cause a significant 
spread in Larmor radii across the data base. In the discharges where ion temperature was 
measured, ion Larmor radius is found to be about 3 mm, no different from many discharges 
with upward sweeping ACs. We then conclude that the downward sweeping ACs are unlikely 
to be due to some anomaly in Larmor radius present in all 26 pulses of the data base. 

Third, we investigate whether the downward sweeping ACs are associated with a specific 
type of plasma equilibrium. With plasma currents varying from 1.1 MA (discharge #68822 
at t = 1.95 s) to 3.5 MA (discharge #46863 at t = 4.8 s), downward sweeping ACs are 
observed with q min « 2, q min « 3, and even q min « 4 (discharge #68822, t = 1.95 s). 
However, equilibrium reconstructions for the discharges where MSE diagnostic data was 
available show that such discharges have very flat g-profiles (See Table [T]) . We note that 
downward sweeping ACs were never observed in JET discharges with strong shear reversal 
and current holes®!. These discharges are typical for JET advanced scenarios and they 
usually exhibit very clear upward sweeping ACs. One can than conjecture that downward 
sweeping ACs require low shear values, which are quite rare in JET discharges with q m in > 1- 



This conjecture is consistent with what follows from Eq. (86) and Fig. |2j shallow g-profiles 
(small q") minimize continuum damping making it easier for the mode to be destabilized 
and grow to an observable level. 

Toroidal mode numbers were determined from Mirnov coil data for 19 discharges con- 
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Pulse # 


tAC S 


Qmin &t tAC 


tMSE S 


(r 2 /g)g" 


Maximum s 


66539 


5.15 


2 


4.96 


0.14 


-0.06 


66550 


5.2 


2 


4.96 


0.4 


-0.1 


72961 


2.7 


3 


2.21 


0.9 


-0.22 


72818 


3.9 


2 


3.06 


0.4 


-0.11 


74896 


2.25 


2 


2.4 


0.35 


-0.07 



TABLE I: Discharges with downward sweeping ACs, for which MSE measurements are available 



taining downward sweeping ACs. All these numbers were positive and lie in the range n = 1 
to n = 5. It is noteworthy that upward sweeping ACs also have only positive toroidal 
mode numbers 6 . The absence of negative toroidal mode numbers agrees with the theoretical 
conclusion that the corresponding modes are stable and damped. 

An idealized model of quiescent plasma with a population of hot ICRF heated ions and a 
shallow (\s\ <C 1) g-profile appears to be relevant to reverse-shear discharges with downward 



sweeping ACs. As follows from Eq. (80), these modes must have positive values of r]. The 
time evolution of r] is dominated by that of nq — m. We thus notionally hold all other 
parameters constant and examine what happens as q (or equivalently 77) decreases with 
time. As the continuum damping decreases with rj the mode eventually becomes unstable. 
If the critical value of rj is above zero we expect to see a downward sweeping trace until rj 
reaches zero whereupon the mode becomes the Geodesic Acoustic Mode 29 . As 77 decreases 
further to negative values, we expect to see an upward sweeping trace until the mode reaches 
the TAE frequency. 

We now consider an example of downward-sweeping ACs shown in Figure [3] Figure [4] 
shows the corresponding power waveforms for the ICRH, NBI and Lower-hybrid current drive 
(LHCD) systems in JET discharge #66550 with toroidal field Bt = 2.46 T and plasma 
current I p i a < 1.9MA It is typical for scenarios aiming at the creation of an internal 
transport barrier that the inductive current drive is still ramping up (see Figure [4]) when the 
main heating power is applied. This creates the non-monotonic g-profile shown in Figure [5] 
(reconstructed from MSE data), with q m in = 2.1 at t = 4.96s when the MSE data was taken. 

This profile has a broad region of very low shear (|s| < 0.1 in the shear-reversed region) 
with r 2 q"/q = 0.4 at the shear reversal point. 
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Pulse No: 66550 Log (IAsint|)l+8) 
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FIG. 3: Spectrogram of O-Mode interferometry from JET pulse #66550 showing ACs (color online) 
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FIG. 4: ICRH, NBI, LHCD power waveforms, and plasma current signal for JET Pulse #66550, 
shaded time interval corresponds to the time shown in Fig. [3] (color online) 



Figure [6] shows the resonant surfaces for the ICRH antennas at the time of observation. 
Most of the RF power was deposited at R = 3.5m. We thus expect a hollow profile of 
fast ions in the region 3m < R < 3.5m. The SELF0 code^ confirms this expectation, 



predicting a hollow profile for the fast particle energy density. As seen from Eq. (73), an 



energetic particle population with such a hollow profile and 9^ < re/ 2 provides a drive for 
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FIG. 5: Reconstructed q profile, from MSE and polarimetry data for JET pulse #66550 
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FIG. 6: ICRH resonant surfaces at time of AC, with reconstructed flux surfaces shown (color 
online) 

AC quasimodes making them observable. 

It is important to emphasize that the quasimodes discussed here do not require a potential 
well, which distinguishes them from the AC eigenmodes discussed in Ref. 05] The drive 
calculated here is destabilizing only for positive toroidal mode numbers, in the typical off- 
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axis ICRH scenarios. Downward sweeping ACs have also been observed in on-axis heated 
discharges. In order to apply our theory to this case we need to include a broader profile 



in bounce angle in our derivation. The generalization of Eq. ( 58 ) to a localized distribution 
with finite width is fairly simple, and one can see quickly that it will be a weighted average 



of Eq. (58) over the range of possible bounce angles. However as K\ is larger in the high field 
side region (by one power of n) we can neglect any contrbution from the low-field-side region 
and conclude that, for on or near axis heating, positive n modes are unstable if J^- < at 
the zero shear point. This agrees with experiment where, unless the zero shear point is very 
close to the magnetic axis, we indeed expect the profile to drop off radially (J^- < 0) at the 
zero shear point. Due to the scarcity of downward sweeping modes on JET, the inability to 
run with hollow profiles and high-field-side IRCH and the scarcity of toroidal mode number 



data it is not possible to completely verify the predictions made after Eq. (73). However 
none of the experimental results on downward sweeping modes from JET contradict our 
predictions. 



VIII. CONCLUSION 



In this paper we have used the action-angle formalism for the bounce-precessional mo- 
tion of the fast ICRH-accelerated ions to calculate the resonant response of the ions to shear 
Alfven perturbations without making the deeply trapped orbit approximation used previ- 
ous jjTiQiii| By exploiting the characteristic features of Alfven Eigenmodes, we reduce the 
complex orbit integrals to local contributions from a few critical points. 

Using this local formulation of the energetic ion response we calculate the resonant ion 
drive for Alfven Cascades (ACs) and AC quasimodes. Based on the expression for the 
drive, and an asymptotic expression for the radiative damping, we conclude that the AC 
quasimode should be observed predominantly in weakly reversed shear configurations and 
that, in such a configuration, it can explain downward sweeping Alfven Cascades below the 
TAE frequency. Furthermore, we conclude that the quasimode should only be excited with 
positive mode numbers in experimentally relevant plasma configurations. 

We reinforce this conclusion via a comparison to experiments carried out on JET. A 
database of JET pulses where quasimodes were observed is compiled and examined, along 
with advanced equilibrium reconstruction to determine the g-profile. The results obtained 
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demonstrate that the observation of downward sweeping Alfven Cascades is indeed linked 
to the presence of a weakly reversed shear plasma. In the few cases where mode number 
information is available only positive mode numbers are excited, agreeing qualitatively with 
our predictions for the resonant excitation mechanism due to trapped energetic ions. 
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APPENDIX A: EVALUATION OF EQS. (49) AND (55) 



In this Appendix we evaluate the integrals K(l,m) and K ext {l } m) defined in Eqs. (49) 



and (55) respectively in the high-n limit. These integrals relate harmonics in straight-field- 
line variables in which the fields are naturally expressed to the action-angle variables (p and 
9 adapted to the particle orbits. 

The major simplifying assumption that allows us to proceed in this appendix is the high 
mode number limit n,m — ► oo, however we immediatly observe that knqR ~ 1 and so 
nq — m doesn't scale with n as n — > oo. We also neglect the small difference 9 — 9 unless it 
is multiplied by the large mode number n or m. 

We note the symmetry in the definitions of K and K ext and so manipulate K ext using 
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Eqs. (p9b and (p| to find 



2n 



dO 



-Ob 



de 
dl 



e -inG{e) X i{f) e -i{nq-m)9 e i T^i u) v-X2{r))e ^_ iW 



Now we can convert this to an integral over 9 

7T 

1 



j^ext _ J_ / jn e -inG(0)xi(r) -i(nq-m)d 1 ^ 
27T 



(Al) 



(A2) 



and comparing this to the definition of in Eq. (49) we find 



K ext (l,m) = K*(l,m) 



with * denoting the complex conjugate. 



(A3) 



Proceeding now with the explicit evaluation of K we use Eq. (49) and Eqs. (29) and (30) 
to obtain 

2ir 7r 



An* 



2tt 



K{l,m) = — / dip / dBeW-^e-*** 



no 



-7T 



e il0 e inG(0)xi(r)+inx2O e i(nq-m 



(A4) 
(A5) 



where I = I — n — . If we are near resonance then / does not scale with n. This integral is 



now of the form: 



x-> 



g(x)e mf{x) dx as n — > oo, 



for which a standard expression exists 
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,(x)c^>dx ~ - 9{X2 ]^ X ^ + ^^if + E 2(c)e^+^)/ 4 . 



2tt 



nf'(x 2 ) nf'(xi) 



c:/'(c)=0 



|n/"(c)|' 
(A6) 



where //(c) = n/"(c)/|n/"(c)|. From this expression we note that such integrals are domi- 
nated by contributions from stationary phase points where f'{x) = and from the ends of 
the interval if no such points exist. We will consider both contributions as stationary phase 
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points do not always exist. The result should be interpreted as providing the leading term. 
If a stationary phase point exists then the endpoint contributions must be ignored. 
To find the stationary phase points, which we denote by 9 = 9*, we solve, 
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which we can rearrange to give 



(G(9)xi(f) + X 29)=0, 



89 xid9 



(A7) 



(A8) 



where t| is related to |j? by Eqs. (36)-(38) , and the ratio Xi/Xu m the low-shear limit, 



becomes 



X2 cos 9b 2PtVL r 



Xi Ro^eQ V M R t 



Eq. (A8) can then be rewritten as 



dG 
~a~9 



cos 9b, 



which using Eq. (28) finally gives 



cos 9* = - cos#b. 
2 



(A9) 



(A10) 



(All) 



If, without loss of generality, we let 9* be the positive solution of this equation we have 
stationary phase points at 9 = ±9*, and thus four stationary phase points in total as each 



9 occurs once on each branch of 9(9) (see Eqs. (30), (34) and (35)). 



Dividing the 9 integral in Eq. (A5) into two intervals, corresponding to the two branches 



of 9(9) (see Eqs. (36) and (37)), we have 



K(l,m) = -(I + J) 



J{nq-m)8 e il8 e inG(e)xi(r)+iriX28 



(A12) 
(A13) 



J 



i(nq-m)6 e %W e inG(6)xi(r)+inx20 



(A14) 



We now calculate I and J in the limit of n — > oo. 



Applying Eq. (A6) to I we find the following contribution from the endpoints 9 = —9b,9 
-7T and 9 = 9 b ,9 = 0: 



e ~i(nq-m)9 b e -iln-inx2K-ixiG(- 



1 c i(nq-m)e b 1 



nX2 



nX2 



(A15) 
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At the stationary phase point, 9*, we have 
d 2 G 



Xi 



d9 2 



^4^?2 f-sinrVcos^-cos^) = G*. (A16) 



e=e , M q 2 R*uj i 



Noting that G* is odd in 9*, we have the contribution to / from 9 = +9* and 9 = —9*, 



e i(nq-m)9* e inG(e*)xi-in(l+ X 2W)+iTTn/4 ^ 



2tt 

(A17) 



\nG 



_|_ e -i(nq-m)6* e inG(-9*)xi-in(l+X2)(ir-6*)-iTT[t/4: _ 



27T 



Where we introduce 9* which is the absolute value of 9 when 9 = 9*, also \i = nG* /\nG*\. 
Thus the two stationary phase points in the interval 9 G [— 7r, 0] are 9 = 9* — n and # = — 9*. 
Combining the above results gives, 

J =c -i(nq-m)0 bc -iin-inx2Tr-inxiG(-0 b ) 1 _ c i(nq-m)6 b 1 (A18) 

nX2 nx2 



+ e 



i(nq-m)e* -i(nx2+l){e*) inG{e*)xi+in^i/A 



27T 



\nG* 

_|_ e -i(nq-m)e* e -i(nx2+i){TT-S*) e inG(-e*)xi-iTrti/4: 



27T 



\nG*\' 

In this expresssion and hereafter G{9) is understood to be evaluated on the branch where 
G(9) < 0. We perform the same analysis for J to give: 

J — e i(nq-m)e b 1 _ c -i(nq~m)6 b c ilTT+inx2^-inxiG{-e b ) 1 (A19) 

nX2 nx2 



+ e 



i(nq-m)6* i{nX2+l) (d*) -inG(6* )xi-infj./4 



27T 



\nG* 

+ e -i(nq-m)9* e i(n X 2+t)(n-e*) e -inG{-6*) X i+iKIJ./± f 



27T 

\nG*[ 



We can then calculate K(l,m) 

K(l,m) = ±-(I + J) 

271 



1 2 

27T U\2 



e - i{nq - m) e b -inxiG(-9 b ) gin ^ + ^^^j 



(A20) 
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To calculate products of K and K* in a simple form we employ an averaging over f, 
which is permitted by Eq. (fi]) as $ doesn't have small scale structure in f. Thus we drop all 
terms containing expressions of the form e m ^ r ^ leaving only those terms which vary slowly 
in the radial direction. This procedure gives, 

K(l,m)K*(l,m + s) = + (e isd * + e~ is6 *) , (A21) 

2ir 2 n 2 xi \niiG*\ x ' 

from which we can calculate K\ and K 2 as defined in 

K 1= Y, m)K™% m + s) = ^ + (A22) 

w 2 n 2 xi n\nG*\ 



=±i 



K 2 =y sK{l, m)K ext (l, m + s) = -i-^^. (A23) 

s=±l Ai 
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